01基因组提取指定区间 - 保留最长200条 - 修改序列名称

从fasta所有序列中提取指定区间

/share/nas1/yuj/script/chloroplast/personality_analysis/dnabarcode/extract_sequence.py

1.从FNA文件中提取指定区间的序列

方法1：使用seqkit工具（推荐）

# 安装seqkit (如果尚未安装)
conda install -c bioconda seqkit

# 提取指定区间 (例如：从100到200的位置)
seqkit subseq --chr "序列ID" -r 100:200 input.fna > output.fna

方法2：使用samtools faidx

# 首先建立索引
samtools faidx input.fna

# 然后提取区间
samtools faidx input.fna "序列ID:100-200" > output.fna

方法3：使用Python脚本

from Bio import SeqIO

def extract_sequence(fna_file, seq_id, start, end, output_file):
    for record in SeqIO.parse(fna_file, "fasta"):
        if record.id == seq_id:
            subseq = record.seq[start-1:end]  # Python是0-based，生物坐标是1-based
            with open(output_file, "w") as f:
                f.write(f">{record.id}_{start}_{end}\n{subseq}\n")
            return
    print(f"序列ID {seq_id} 未找到")

# 使用示例
extract_sequence("input.fna", "chr1", 100, 200, "output.fna")

方法4：使用bedtools

# 首先创建一个BED文件定义区间
echo -e "序列ID\t99\t200" > region.bed  # BED是0-based半开区间

# 然后提取序列
bedtools getfasta -fi input.fna -bed region.bed -fo output.fna

注意事项

  1. 区间定义：不同工具对区间定义可能不同(0-based或1-based，闭区间或半开区间)
  2. 序列ID：确保输入的序列ID与文件中的完全匹配
  3. 大文件处理：对于大基因组文件，建议先建立索引(如samtools faidx)
  4. 反向互补链：如果需要提取反向互补链，可以添加参数如--reverse-complement(seqkit)或处理Python中的序列

2.从FASTA文件 map_gene.fa 中保留最长的200条序列

方法2：使用 seqkit（更高效，推荐）

seqkit sort --by-length --reverse map_gene.fa \
  | seqkit head -n 200 > top200.fa

步骤解释：

  1. seqkit sort 按序列长度降序排序。
  2. seqkit head 提取前200条记录。

安装 seqkit：

wget https://github.com/shenwei356/seqkit/releases/download/v2.3.0/seqkit_linux_amd64.tar.gz
tar -zxvf seqkit_linux_amd64.tar.gz
sudo mv seqkit /usr/local/bin/

方法3：纯 awk + 基础命令

awk '/^>/ { 
        if (header) print header ORS seq
        header=$0
        seq=""
        next
     } 
     { seq = seq $0 } 
     END { if (header) print header ORS seq }' map_gene.fa \
  | awk 'BEGIN {RS=">"; FS="\n"} NR>1 { 
        len=0; seq=""
        for (i=2; i<=NF; i++) { len+=length($i); seq=seq $i }
        print len, ">"$1, seq
     }' \
  | sort -k1,1nr \
  | head -n 200 \
  | cut -d' ' -f2- \
  | tr ' ' '\n' \
  | awk '/^>/ {print; next} {gsub(/.{80}/,"&\n"); printf "%s", $0}' > top200.fa

步骤解释：

  1. 合并多行序列为单行。
  2. 计算每条序列的长度并附加到行首。
  3. 按长度降序排序，取前200条。
  4. 恢复FASTA格式，每80字符换行。

3.批量修改FASTA文件序列ID

方法1：使用Linux命令 awk

awk 'BEGIN {
    # 加载chrom_map.txt到关联数组map中
    while (getline < "chrom_map.txt") {
        map[$1] = $2
    }
}
# 处理以">"开头的行（序列ID行）
/^>/ {
    # 提取旧ID（去掉开头的">"，取第一个字段）
    old_id = substr($1, 2)
    # 如果旧ID在映射表中，则替换为新ID
    if (old_id in map) {
        $1 = ">" map[old_id]  # 修改第一个字段为新ID
    }
}
# 打印所有行（1表示默认动作：打印当前行）
1' input.fasta > output.fasta